C.................... RSPE2B.FOR ...............................
C...... A program to calculate photoelectron fluxes- by Phil Richards
C...... 2-stream interhemispheric photoelectron program loosely based on
C...... a program by Nagy and Banks but incorporating several of my own
C...... innovations. see publications. a more accurate model is available
C...... on <pe2s.for> but this one is accurate enough for most purposes.
C---      Z= altitude, ZOX = [O], ZN2= [N2], ZO2= [O2], HE= [He]
C---      BM= magnetic field strength, BM= field l,predine coordinate factor
C---      IMAX= # of grid points on field line, DH= distance between
C---      points in modified (p,q) system. N= density of H+, O+, minor
C---      ions(O2+ + NO+), and He+. TI= temperature of H+, O+, electrons
C---      EHT(3,I)= electron heating rate to be returned
C---      TN= neutral temperature, SZA= solar zenith angle, 
C---      FRPAS= fraction of flux lost in plasmasphere
C---      IVLPS= solution switch, not important. ZWR= altitude for
C---      printing spectrum. ISOL= solution switch, not important
C---      Double Precision could be switched off
      SUBROUTINE PE2S(Z,ZOX,ZN2,ZO2,HE,BM,BG,IMAX,DH,N,TI,EHT
     > ,TN,SZA,FRPAS,IVLPS,ZWR,ISOL,SL)
      DOUBLE PRECISION Z,ZOX,ZN2,ZO2,BM,BG,DH,N,TI,EHT,TN,SZA,FRPAS
     >  ,HE(401),FD(9),TEFLUX,TFLUX,AVEE,TEEHT,FACPAS,GL,SL(401)
     >  ,COLUM,OTHPR1,OTHPR2,OTHPR3
      REAL AP,DEC,ETRAN,BLON,F107,F107A,SEC,ZN2D
      REAL EFAC,ECDT,ECHROM,ECORON    !.. Solar eclipse factors and time step
      REAL ECFACN(901),ECFACS(901)    !.. Eclipse factors for photoelectrons
      INTEGER FEXISTS,IWR,IMAX1,IMXM,IPASC,IEQ,IPAS,ICON,IWMIN
      DIMENSION TSIGNE(401),TN(401),BG(401),SECSAV(2,401),IE(201)
     > ,SIGEX(3),SIGION(3),E(201),SZA(401),Z(401),XN(3,401),FYSUM(401)
     > ,EHT(3,401),SPRD(3,6),ZOX(401),ZN2(401),ZO2(401),BM(401)
     > ,N(4,401),TI(3,401),PRED(401),PRODWN(201,401),EQN2D(401)
     > ,PRODUP(201,401),EB(201),ALTA(401),DELTE(201),JOE(201)
     > ,JN2E(201),PROB(201),PARSIG(22),SIGEL(3),PEBSC(3)
      REAL PRIMION(401),PETION(401),HE10830XS
      COMMON/TRI/PHIDWN(401),PHIUP(401),T1(401),T2(401),DS(401)
      COMMON/PANGS/YPAS,PASK,IPAS,IPASC,FPAS
      !-- Common to hold the EUV and photoelectron production rates
      COMMON/EUVPRD/EUVION(3,12,401),PEXCIT(3,12,401),PEPION(3,12,401)
      COMMON/EFLUX/TEFLUX,TFLUX,AVEE,TEEHT,FACPAS
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/PE2COM/PROB    !.. common suggested by Kent Miller
      COMMON/NPRD/COLUM(3,401),OTHPR1(6,401),OTHPR2(6,401),OTHPR3(6,401)
      COMMON/SOL/UVFAC(59),EUV

      DATA  M,ZLB,ZPROD,JMIN,M1,AVMU,EMIN,EMAX,ECHGE,IU9,ZPAS
     >  / 41, 120., 999., 2  ,1 ,.577,1.0 , 800., 30. , 9 ,1000./
      DATA ICON /0/,JMAX/0/,EBSCX/0.0/
      DATA M2/0/,IPR/0/,IWR/0/,JN2E/201*0/,SIGEX/3*0.0D0/
      DATA DELTE/201*0.0/
      DATA JOE/201*0/,SIGION/3*0.0D0/,EB/201*0.0/,E/201*0.0/
      DATA XN/1203*0.0/, ZN2D/180./

      !.. turns on different flux prints, 0= up+down fluxes, 1= up and 
      !.. down fluxes, 2= cross sections etc., 3=up and down for FAST
      DATA IPFLUX/1/

      DATA EPOT,ELOSSOX,ELOSSN2,ELOSSO2,SHAPE/17.0,11.0,10.0,9.0,14.0/
      !.. Burnett&Rountree O branching ratios. O2. need revision
      !.. Shemansky and Liu N2 cross sections, JGR 2005
      DATA SPRD/.4,.56,.44, .4,.28,.44, .2,.06,.10, 0.,.05,.00, 0.,.05
     > ,.00, 0.0,0.0,0.02/
      DATA IE/201*0/, F107SV/0.0/,FEXISTS/-9/

      !.. Setting up energy cell boundaries first time thru only
      !.. If Printing, go to 1 eV resolution  
      IF(JMAX.LE.0)  THEN
         IF(ZWR.LT.80) THEN
           !.. can use lower resolution for the FLIP run
           CALL ECELLS(35.0,800.0,1.0,50.0,EMAX,JMAX,EB,E,DELTE)
         ELSE
           !.. use high resolution for printing
           CALL ECELLS(65.0,800.0,1.0,20.0,EMAX,JMAX,EB,E,DELTE)
         ENDIF
         !.. set up bins for degraded primaries from ionizations
         DO J=JMAX,2,-1
           !.. Calculate the average energy of the secondaries.
           ELIM=0.5*(E(J)-EPOT)  !.. Maximum secondary energy
           !.. normalize secondary distribution
           FNORM=1.0/ATAN(ELIM/SHAPE)/SHAPE     
           AVESEC=FNORM*0.5*SHAPE**2*ALOG(1.0+(ELIM/SHAPE)**2)
           ELOSS=EPOT+AVESEC

           !.. allocate bins for degraded primaries for ionization, and
           !.. excitation of O and N2
           IF(E(J).GT.ELOSS) CALL FNDBIN(J,JMAX,ELOSS,E,IE)
           IF(E(J).GT.ELOSSOX) CALL FNDBIN(J,JMAX,ELOSSOX,E,JOE)
           IF(E(J).GT.ELOSSN2) CALL FNDBIN(J,JMAX,ELOSSN2,E,JN2E)
         ENDDO

         !.. proportion of total secondary ions in energy bin E(J). Note that 
         !.. a running sum of secondary electrons is kept but only distributed
         !.. immediately prior to calculating the flux at the next lowest energy. 
         !.. The apportionment is based on an empirical calculation using the 
         !.. full Opal et al. [1971] distributions
         DO J=1,JMAX
           PROB(J)=0.2526*EXP(-0.2526*E(J))
         ENDDO
      ENDIF

      !........ Set up production frequencies for the energy cells
      !.. see if EUV data file exists and then calculate frequencies
      IF(FEXISTS.EQ.-9) CALL TFILE(42,FEXISTS)
      IF(FEXISTS.NE.0)
     >  CALL CONV_FILE_FREQ(JMAX,F107,F107A,SEC/3600.0,E,DELTE)

      !... Use model fluxes
      IF(FEXISTS.EQ.0.AND.ABS((F107-F107SV)/F107).GE.0.05) THEN
        F107SV=F107

        IF(NINT(UVFAC(58)).EQ.-1.OR.UVFAC(58).EQ.-3) THEN
            CALL CONVF_EUVAC(JMAX,F107,F107A,E,DELTE)	!.. Standard EUVAC F107 scaling
        ELSEIF(NINT(UVFAC(58)).EQ.-2.OR.UVFAC(58).EQ.-4) THEN
            CALL CONVF_SEE_V8(JMAX,F107,F107A,E,DELTE)  !.. SEE F107 scaling
           !.. The following are for testing various EUV fluxes 
           !..CALL CONV_SC21REF(JMAX,F107,F107A,E,DELTE) !.. HFG F107 scaling
           !..CONV_FILE_FREQ is used for testing one-off spectrum. No F107 scaling
         ENDIF !.. Endif TFILE
      ENDIF  !.. Endif F107

      TEEHT=0.0
      TEFLUX=0.0
      TFLUX=0.0
      IEQ=(IMAX+1)/2

      !.. Set pitch angle trapping factor. FPAS is controlled by input
      !.. parameter FRPAS. There are several options. NOTE - to find
      !.. code related to pitch angle scattering search for string PAS
      FPAS=0.0
      !.. Standard - user specifies 0.0 <= FPAS <= 1.0 
      IF(FRPAS.GE.0.0.AND.FRPAS.LE.1.0) THEN
        FPAS=FRPAS
      !.. If FRPAS<0, adjust fraction of trapping according to TEC
      ELSEIF(FRPAS.LT.0) THEN
        FPAS=DMIN1(DABS(FRPAS),6.0E-6*DABS(FRPAS)*N(2,IEQ)*
     >     ((1+Z(IEQ)/6370.)**4-(1+ZPAS/6370.)**4))
      !.. If FRPAS > 1 Stop electrons in plasmasphere but no extra heating
      ELSEIF(FRPAS.GT.1.0) THEN
        FPAS=1.0
      ENDIF
      !.. Make sure 0.0 <= FPAS <= 1.0 
      IF(FPAS.GT.1.0) FPAS=1.0
      IF(FPAS.LT.0.0.OR.Z(IEQ).LE.ZPAS) FPAS=0.0

      DO 5 I=1,IMAX
      DO 5 IK=1,12
      DO 5 IS=1,3
      PEPION(IS,IK,I)=0.0
      PEXCIT(IS,IK,I)=0.0
 5    CONTINUE

      DO 6 I=M1,IMAX
      DO 6 J=1,JMAX
      PRODUP(J,I)=0.0
 6    PRODWN(J,I)=0.0

      DO I=M1,IMAX
      SECSAV(1,I)=0.0
      SECSAV(2,I)=0.0
      PHIUP(I)=0
      PHIDWN(I)=0
      ALTA(I)=Z(I)
      XN(1,I)=ZOX(I)
      XN(2,I)=ZO2(I)
      XN(3,I)=ZN2(I)
      EHT(3,I)=0.0D0
      PETION(I)=0.0
      ENDDO

      !.. trap large printing altitude
      IF(ZWR.GT.10*6370.00.OR.ZWR.GT.Z(IMAX/2+1)) THEN
      ZWR=Z(IMAX/2)
      WRITE(IU9,'(/A,1P,E12.2/)') 
     >    ' *** Photoelectron print altitude too large, reset to =',ZWR
        WRITE(8,'(/A,1P,E12.2/)') 
     >    ' *** Photoelectron print altitude too large, reset to =',ZWR
      ENDIF
      
      !........ Calculate arc length DS
      DO 13 I=M1,IMAX
      IF(I.GT.1) DS(I)=SL(I)-SL(I-1)
      IF(I.GT.1) DS1=0.5*(DH/BG(I)+DH/BG(I-1))
      IF(Z(I).LT.ZPAS.AND.I.LT.IEQ) IPAS=I
      IF(ZWR.GT.80.AND.Z(I).LT.ZWR.AND.I.LE.IEQ) IWR=I
      IF(Z(I).GT.ZPROD) GO TO 13
      IF(Z(I).LT.400.AND.I.LT.IEQ)  M400=I
      IF(I.LT.IEQ)  M=I
      IF(Z(I).LT.ZLB.AND.I.LT.IEQ)  M2=I
13    CONTINUE

      !..WRITE(6,'(2I,2F7.1,1P,9E10.2)') M1,IWR,Z(M1),Z(IWR),BM(M),BM(IWR)
      IF(ZLB.LE.Z(2)) M2=2
      IMAX1=2*IEQ-M2
      IMXM=IMAX+1-M
      !.. make sure closest altitude to ZWR is selected for printing
      IF(ZWR.GT.80) THEN
        IF(ABS(Z(IWR+1)-ZWR).LT.ABS(ZWR-Z(IWR))) IWR=IWR+1
      ENDIF
      ICON=2*IEQ-IWR     !.. alt. index of conjugate to IWR
      IWMIN=MAX(IWR-6,1) !.. alt. index of lowest altitude to print
      IPASC=2*IEQ-IPAS   !.. alt. index for pitch angle scattering
     
      !....... write altitudes for photoelectron fluxes
      IF(ZWR.GT.80.AND.ZWR.GE.Z(1)) THEN
         !.. For up and down fluxes
         IF(IPFLUX.EQ.1) THEN
           WRITE(IU9,316) (NINT(Z(K)),K=IWMIN,IWMIN+9),
     >       (NINT(Z(K)),K=IWMIN,IWMIN+9)
           WRITE(8,316) (NINT(Z(K)),K=IWMIN,IWMIN+9),
     >       (NINT(Z(K)),K=IWMIN,IWMIN+9)
         ELSEIF(IPFLUX.EQ.3) THEN
           WRITE(IU9,317) (NINT(Z(K)),K=IWMIN,IWMIN+9),
     >       (NINT(Z(K)),K=IWMIN,IWMIN+9)
           WRITE(8,317) (NINT(Z(K)),K=IWMIN,IWMIN+9),
     >       (NINT(Z(K)),K=IWMIN,IWMIN+9)
         !.. for total fluxes
         ELSE IF(IPFLUX.EQ.0) THEN
           WRITE(IU9,318) (NINT(Z(K)),K=IWMIN,IWMIN+9)
           WRITE(8,318)   (NINT(Z(2*IEQ-K)),K=IWMIN,IWMIN+9)
         ENDIF
      ENDIF

 316   FORMAT(/'  *** Photoelectron fluxes (eV cm2 str s1) for'
     >   ,1X,'several altitudes (km)',/47X,
     >   '<--UPWARD-->',36X,'-||-',35X,'<--DOWNWARD-->',/3X,'E',33I9)

 317   FORMAT(/'  *** Photoelectron fluxes (eV cm2 str s1) for'
     >   ,1X,'several altitudes (km) scaled for comparison with FAST'
     >   ,/47X,
     >   '<--UPWARD-->',36X,'-||-',35X,'<--DOWNWARD-->',/3X,'E',33I9)

 318   FORMAT(/'    Total photoelectron fluxes (phiup+phidown)
     > /(eV cm2 str s1) for a number of altitudes (km)',/3X,'E',33I9)

      IF(IPFLUX.EQ.2) THEN
        WRITE(IU9,'(/5X,A,F10.2,A)') 
     >   ' e* fluxes, X-sects, etc. at alt =',Z(IWR),' km'
        WRITE(8,'(/5X,A,F10.2,A)')
     >   ' e* fluxes, X-sects, etc. at alt =',Z(ICON),' km'
      ENDIF
      IF(IPFLUX.EQ.2) WRITE(IU9,319)
      IF(IPFLUX.EQ.2) WRITE(8,319)
 319  FORMAT('  E      SIGEXO  SIGIONO  SIGEXN2  SIGION2  SIGEL'
     > '    PRED     FYSUM    TSIGNE   PHIUP    PHIDWN   PRODUP'
     > '   PRODWN  SEC_UP   SEC_DWN')

      !......... Determine flux tube volume for heating due to pitch angle
      !......... trapping
      VTOT=0.0
      DO 413 I=IPAS,IPASC-1
 413  VTOT=VTOT+DS(I)/(BM(I)*2.038E+8)
      IF(VTOT*BM(M2).GT.0.0) PASK=FPAS/(VTOT*2.038E+8*BM(M2))

      !.. Calculate solar ECLIPSE factors if any. See EFAC near the PEPRIM call
      CALL ECFUN(SEC,IDAY,45.0,ECDT,ECHROM,ECORON)
      CALL PEFACS(ECORON,ECHROM,ECFACN) !.. adjust the PE factors for 304A        
      CALL ECFUN(SEC,IDAY,-45.0,ECDT,ECHROM,ECORON)
      CALL PEFACS(ECORON,ECHROM,ECFACS) !.. adjust the PE factors for 304A        

      J=JMAX      !.......... highest energy set for J

C////////////main calculations  begin here ////////////
 23   CONTINUE

      !-- Pitch angle trapping from Khazanov et al. 1992
      IF(FRPAS.LT.0) THEN
          FPAS=-FRPAS*(-0.016667*E(J)+0.9167)
          IF(FPAS.GT.1.0) FPAS=1.0
          IF(FPAS.LT.0.0) FPAS=0.0
          IF(VTOT*BM(M2).GT.0.0) PASK=FPAS/(VTOT*2.038E+8*BM(M2))
      ENDIF

      !.......... Evaluate ionization branching ratios for O+
      CALL OXRAT(E(J),SPRD(1,1),SPRD(1,2),SPRD(1,3))

      !...  O, O2, and N2 elastic cross sections from AUR2SA MAY 1999
      CALL ELASTC(E(J),PEBSC,SIGEL)

      !....... set up the probability distribution for the secondary electrons
      !....... T1, T2 are coeefs of d.e. b&k p258 ............
      DO 36 I=M1,IMAX
      PRED(I)=0.0D0
36    CONTINUE

      !............. Loop for calculating primary and cascade production
      DO 37 I=M1,IMAX
      IF(Z(I).GT.ZPROD) GO TO 35
      IF(I.LT.JMAX/2) EFAC=ECFACN(J)   !.. set north eclipse function
      IF(I.GE.JMAX/2) EFAC=ECFACS(J)   !.. set south eclipse function
      JMONO=-177           !.. index for mono-energetic fluxes
      IF(J.ne.JMONO)
     > CALL PEPRIM(I,XN,HE(I),PRED,JMAX,J,E(J),DELTE(J),IWR,EFAC)
      !......... add electron quenching of N(2D) and O+(2D) to primary prod
        PRIMION(I)=PRIMION(I)+PRED(I)*DELTE(J)  !.. primary ion production
        IF(Z(I).GT.ZN2D.AND.ISOL.NE.0.AND.E(J).LE.4.AND.E(J).GT.3) THEN
           CALL CMINOR(0,I,0,IPR,FD,11,HE(I))
           EQN2D(I)=FD(8)
           PRED(I)=PRED(I)+FD(6)
        ENDIF
        IF(Z(I).GT.ZN2D.AND.ISOL.NE.0.AND.E(J).LE.3.AND.E(J).GT.2)
     >       PRED(I)=PRED(I)+EQN2D(I)
        IF(I.EQ.-IWR) WRITE(6,'(2F9.2,1P,2E10.2,0P,F9.2,2I7)') 
     >  Z(I),E(J),PRED(I),EQN2D(I)
     > ,DELTE(J),INT(1+E(J)-0.5*DELTE(J)),INT(E(J)+0.5*DELTE(J))
      !...... Total energy deposition to photoelectrons
      EUVION(1,11,I)=EUVION(1,11,I)+PRED(I)*E(J)*DELTE(J)
      !....... cascade production using frequencies ....
 35   CONTINUE
      IF(Z(I).GT.ZPROD) GO TO 37
      !....     IF(E(J).LE.3.AND.E(J).GT.2)
      !....    >  WRITE(3,313) Z(I),PRED(I),FD(8),PRODWN(J,I),PRODWN(J,I)
 37   CONTINUE

      !.. Get total cross sections
      CALL SIGEXS(E(J),SIGEX,SIGION,SIGO1D) !.. PGR cross sections

      !.. CALL STRICKXS(J,JMAX,E,SIGION,SIGEX,SIGO1D)  !.. Strickland cross sections
      !..WRITE(6,'(I5,F10.2,1P,9E10.2)') J,E(J),SIGION(1),SIGEX(1)
      !..>   ,SIGION(3),SIGEX(3),SIGO1D

      !.. GLOW total cross sections
      !..CALL GLOWSIGS(E(J),SIGEX,SIGION)

      !.. Get OX partial cross sections
      CALL OXSIGS(E(J),PARSIG,TSIG)

      !....... energy loss to thermal electrons by coulomb collisions
      DO 29 I=M1,IMAX
        ET=8.618E-5*TI(3,I)
        ZNE=N(1,I)+N(2,I)+N(3,I)+N(4,I)
        TSIGNE(I)=((3.37E-12*ZNE**0.97)/(E(J)**0.94))
     > *((E(J)-ET)/(E(J)-(0.53*ET)))**2.36/(E(J)-E(J-1))
        T1(I)=0.0
29    CONTINUE

      !........ See if energy loss exceeds bin size
      SUMTSE=0.0
      DE=0.0
      DO 11 I=M1,IMAX
        TEST=TSIGNE(I)*DS(I)
        IF(Z(I).GT.ZPROD) SUMTSE=SUMTSE+TEST
        IF(Z(I).GT.Z(IPAS)) DE=DE+TEST
        IF(TEST.LE.1.0) GO TO 11
          TSIGNE(I)=1.0/DELTE(J)/DS(I)
 11   CONTINUE

      !....... t1,t2 are coeefs of d.e. b&k p258 ............
      !....... Calculate total scattering cross section
      DO 30 I=M1,IMAX
      T2(I)=TSIGNE(I)
      DO 30 IS=1,3
      T1(I)=T1(I)+XN(IS,I)*SIGEL(IS)*PEBSC(IS)
      TABSXS=SIGEX(IS)+SIGION(IS)
      T2(I)=T2(I)+XN(IS,I)*(SIGEL(IS)*PEBSC(IS)+TABSXS)
 30    CONTINUE


      !.......   fluxes in local equil and set boundary conditions on fluxes
      DO 54 ICN=1,M2
        ICF=2*IEQ-ICN
        PHIDWN(ICN)=(.5*PRED(ICN)+PRODWN(J,ICN))/(T2(ICN)-T1(ICN))
        PHIDWN(ICF)=(.5*PRED(ICF)+PRODWN(J,ICF))/(T2(ICF)-T1(ICF))
        PHIUP(ICN)=PHIDWN(ICN)
        PHIUP(ICF)=PHIDWN(ICF)
 54   CONTINUE

      !........... take the adiabatic variation of pitch angle into account
      CALL PITCH(J,0,M,M1,IMAX,T1,T2,PRED,PRODUP,PRODWN,BM,Z)

      !....... calculate interhemispheric fluxes ......
      DO 65 ITS=1,4
      CALL TRIS1(1,M2,M,J,M,BM,BG,Z,IMAX,DH,PRED,PRODUP,PRODWN)
      CALL TRISM1(-1,IMXM,IMAX1,J,M,BM,BG,Z,IMAX,DH,PRED,PRODUP,PRODWN)
 65    CONTINUE
      !.......... reset the fluxes
      CALL PITCH(J,1,M,M1,IMAX,T1,T2,PRED,PRODUP,PRODWN,BM,Z)

      !...... EHPAS=heating due to pitch angle trapping: DE= normal loss in the
      !...... protonosphere to electrons. Modification made 1/29/1996 ->
      !...... EHPAS= Energy * (FPAS/VOL) * flux * delta E * area at PAS alt.

      EHPAS=E(J)*PASK*(PHIUP(IPAS)+PHIDWN(IPASC))*DELTE(J)*
     >       BM(M2)/BM(IPAS)
      IF(DE.GT.E(J).OR.NINT(FRPAS).EQ.2) EHPAS=0.0
      IF(BM(IPAS).GT.0.0) TEFLUX=
     > TEFLUX+E(J)*(PHIUP(IPAS)+PHIDWN(IPASC))*DELTE(J)*BM(M2)/BM(IPAS)
      IF(BM(IPAS).GT.0.0)
     > TFLUX=TFLUX+(PHIUP(IPAS)+PHIDWN(IPASC))*DELTE(J)*BM(M2)/BM(IPAS)

      !...... cross section for O(1S) Jackman et al. 1977
      SIGO1S=0.0
      IF(E(J).GT.4.17) SIGO1S=6.54E-17*(1-SQRT(4.17/E(J)))/E(J)
      !......... Borst and Zipf 3914 cross section used to calculate the branching 
      !......... ratio to the B state. Factor 1.54 is inverse Frank-Condon
      S3914=0.0
      SPRD(3,3)=0.0
      IF(E(J).GT.19.0) THEN
        S3914=8.83E-16*(1.0-9.0/E(J))**7.47/E(J)**0.7
        SPRD(3,3) = 1.54 * S3914/SIGION(3)
      ENDIF

      !........ Save total flux and fluxes for cascade calculation
      DO 19 I=M1,IMAX
      FYSUM(I)=(PHIUP(I)+PHIDWN(I))
 19   CONTINUE

      !.. Calculate thermal electron heating and ion production rates
      DO 20 I=M1,IMAX
      PHISUM=FYSUM(I)*DELTE(J)
      !.. electron heating. Add extra heat from pitch angle trapping
      EHT(3,I)=EHT(3,I)+PHISUM*TSIGNE(I)
      IF(Z(I).GT.Z(IPAS)) EHT(3,I)=EHT(3,I)+EHPAS
      IF(Z(I).GT.Z(M)) GO TO 20

      !.. excitation of O, O2, and N2
      DO 21 K=1,6
 21   PEXCIT(1,K,I)=PEXCIT(1,K,I)+PARSIG(K)*(PHISUM*XN(1,I))

      PEXCIT(1,12,I)=PEXCIT(1,12,I)+(SIGEX(1)-SIGO1D)*(PHISUM*XN(1,I))
      PEXCIT(2,12,I)=PEXCIT(2,12,I)+SIGEX(2)*PHISUM*XN(2,I)
      OTHPR2(4,I)=OTHPR2(4,I)+HE10830XS(E(J))*PHISUM*HE(I)  !.. He 10830 A prod
      !..OTHPR1(2,I)=OTHPR1(2,I)+0.1*SIGION(3)*PHISUM*HE(I) !.. He+ from PE
      CALL EPN2XS(I,XN,PHISUM,ALTA(I),E(J),SIGEX,SIGION,PEXCIT)

      !.. secondary ionization.
 22   DO 130 IS=1,3
      PRION=PHISUM*SIGION(IS)*XN(IS,I)
      PETION(I)=PETION(I)+PRION
      DO 25 IK=1,6
         PEPION(IS,IK,I)=PEPION(IS,IK,I)+PRION*SPRD(IS,IK)
 25   CONTINUE
 130  CONTINUE
 20   CONTINUE

      !.. printing fluxes
      IF(ZWR.GT.80) THEN
        !.. Components of equation
        IF(IPFLUX.EQ.2) THEN
          WRITE(IU9,313) E(J),SIGEX(1),SIGION(1),SIGEX(3),SIGION(3)
     >    ,SIGEL(1),PRED(IWR),FYSUM(IWR),TSIGNE(IWR),PHIUP(IWR)
     >    ,PHIDWN(IWR),PRODUP(J,IWR),PRODWN(J,IWR)
     >    ,SECSAV(1,IWR),SECSAV(2,IWR)
          WRITE(8,313) E(J),SIGEX(1),SIGION(1),SIGEX(3),SIGION(3)
     >    ,SIGEL(1),PRED(ICON),FYSUM(ICON),TSIGNE(ICON),PHIUP(ICON)
     >    ,PHIDWN(ICON),PRODUP(J,ICON),PRODWN(J,ICON)
     >    ,SECSAV(1,ICON),SECSAV(2,ICON)        !.. Up and down fluxes
        ELSEIF(IPFLUX.EQ.1) THEN
          WRITE(IU9,313) E(J),(PHIUP(K)/6.3,K=IWMIN,IWMIN+9),
     >      (PHIDWN(K)/6.3,K=IWMIN,IWMIN+9)
          WRITE(8,313) E(J),(PHIDWN(2*IEQ-K)/6.3,K=IWMIN,IWMIN+9)
     >      ,(PHIUP(2*IEQ-K)/6.3,K=IWMIN,IWMIN+9)
        !.. Up and down fluxes for comparison with Peterson data. Note that 
        !.. the magnetic field factor has to be removed to agree with data
        !.. The scaling is done from 400 km to the satellite altitude.
        ELSEIF(IPFLUX.EQ.3) THEN
          WRITE(IU9,313) E(J)
     >      ,(PHIUP(K)*BM(M400)/BM(K)/6.3,K=IWMIN,IWMIN+9)
     >      ,(PHIDWN(K)*BM(M400)/BM(K)/6.3,K=IWMIN,IWMIN+9)
          WRITE(8,313) E(J)
     >      ,(PHIDWN(2*IEQ-K)*BM(M400)/BM(2*IEQ-K)/6.3,K=IWMIN,IWMIN+9)
     >      ,(PHIUP(2*IEQ-K)*BM(M400)/BM(2*IEQ-K)/6.3,K=IWMIN,IWMIN+9)
        !.. total fluxes
        ELSEIF(IPFLUX.EQ.0) THEN
          WRITE(IU9,313) E(J),
     >      ((PHIUP(K)+PHIDWN(K))/12.57,K=IWMIN,IWMIN+9)
          WRITE(8,313) E(J),
     >      ((PHIDWN(2*IEQ-K)+PHIUP(2*IEQ-K))/12.57,K=IWMIN,IWMIN+9)
        ENDIF
      ENDIF

      !.......... Calculate Cascade production ....
      DO 83 I=M1,IMAX
      !.. Cascade from thermal electron collisions
      TSIGNE(I)=TSIGNE(I)*DELTE(J)/DELTE(J-1)
      PRODUP(J-1,I)=PRODUP(J-1,I)+PHIUP(I)*TSIGNE(I)*(1-EBSCX)
     > +PHIDWN(I)*TSIGNE(I)*EBSCX
      PRODWN(J-1,I)=PRODWN(J-1,I)+PHIDWN(I)*TSIGNE(I)*(1-EBSCX)
     > +PHIUP(I)*TSIGNE(I)*EBSCX

      IF(Z(I).GT.ZPROD) GO TO 83

      !.. calculate secondary and cascade production from ionizing
      !.. collisions
      IF(J.GT.0) THEN
        CALL CASION(J,I,IWR,ALT,XN,PRODUP,PHIUP(I),PHIDWN(I)
     >   ,E,DELTE,1,SECSAV,SIGION,JMAX,IE(J),ELOSS,PROB(J))
        CALL CASION(J,I,IWR,ALT,XN,PRODWN,PHIDWN(I),PHIUP(I)
     >   ,E,DELTE,2,SECSAV,SIGION,JMAX,IE(J),ELOSS,PROB(J))
      ENDIF
      !.. calculate cascade from atomic O - O(1D)
      CALL CASEX(J,1,I,IWR,ALTA(I),SIGEX,XN,PRODUP,PHIUP,PHIDWN
     >  ,E,DELTE,SIGO1D,SIGION,JMAX,JOE(J),ELOSSOX)
      CALL CASEX(J,1,I,IWR,ALTA(I),SIGEX,XN,PRODWN,PHIDWN,PHIUP
     >  ,E,DELTE,SIGO1D,SIGION,JMAX,JOE(J),ELOSSOX)

      !.. cascade from N2
      CALL CASEX(J,2,I,IWR,ALTA(I),SIGEX,XN,PRODUP,PHIUP,PHIDWN
     >   ,E,DELTE,SIGO1D,SIGION,JMAX,JN2E(J),ELOSSN2)
      CALL CASEX(J,2,I,IWR,ALTA(I),SIGEX,XN,PRODWN,PHIDWN,PHIUP
     >   ,E,DELTE,SIGO1D,SIGION,JMAX,JN2E(J),ELOSSN2)

      !.. cascade from O2
      CALL CASEX(J,3,I,IWR,ALTA(I),SIGEX,XN,PRODUP,PHIUP,PHIDWN
     >   ,E,DELTE,SIGO1D,SIGION,JMAX,JN2E(J),ELOSSN2)
      CALL CASEX(J,3,I,IWR,ALTA(I),SIGEX,XN,PRODWN,PHIDWN,PHIUP
     >   ,E,DELTE,SIGO1D,SIGION,JMAX,JN2E(J),ELOSSN2)

      !... Cascade from O(1D) and N2(v)
      CALL CASSIM(J,I,IWR,ALTA(I),SIGEX,XN,PRODUP,PHIUP,PHIDWN
     >   ,E,DELTE,SIGO1D,SIGION,JMAX)
      CALL CASSIM(J,I,IWR,ALTA(I),SIGEX,XN,PRODWN,PHIDWN,PHIUP
     >  ,E,DELTE,SIGO1D,SIGION,JMAX)
 83   CONTINUE
 82   J=J-1

      IF(E(J).LT.EMIN) GO TO 80
      IF (J .LT. 2) GO TO 80
      GO TO 23
      !............................. END OF MAIN ENERGY LOOP ..............

 80   CONTINUE

      !... add in energy from last energy bin
      DO 180 I=1,IMAX
         ELEFT=(PRODUP(1,I)+PRODWN(1,I))*E(1)
         EHT(3,I)=EHT(3,I)+ELEFT
         IF(Z(I).GT.ZN2D.AND.Z(I).LT.ZPROD) THEN
             CALL CMINOR(0,I,0,IPR,FD,12,HE(I))
             EPHOT=EHT(3,I)-FD(8)*2.4-ELEFT
         ENDIF
         RNENI=(N(1,I)+N(2,I)+N(3,I))*(N(1,I)/16.0+N(2,I)+N(3,I)/30.0)
         HKEI=7.7E-6*RNENI*(TI(3,I)-TI(2,I))/(TI(3,I)**1.5)
         IF(Z(I).GT.Z(IPAS)) TEEHT=TEEHT+
     >     (EHT(3,I)-HKEI)*DS(I)*BM(IPAS)/BM(I)
 180  CONTINUE
      IF(J.NE.-999) RETURN  
      IF(TFLUX.GT.0) AVEE=TEFLUX/TFLUX
      !.. This section for calculating the energy loss per ionization
      IF(ZWR.GT.80) THEN
      IF(JMONO.GT.0)  WRITE(IU9,'(A,F8.2)') 'Mono Energy=', E(JMONO)
        WRITE(IU9,'(A,A)') '  ALT    EV/ION  PE_ION EUV_ION TOT_ION'
     >     ,'    EUV/PE    Etot'
        DO I=2,IWR
          WRITE(IU9,'(5F8.1,F10.3,1P,9E10.2)') Z(I)
     >      ,EUVION(1,11,I)/PETION(I),PETION(I),PRIMION(I)
     >      ,PETION(I)+PRIMION(I),PRIMION(I)/PETION(I),EUVION(1,11,I)
        ENDDO
      ENDIF
      RETURN

 95   FORMAT(2X,I6,F8.1,1P,22E12.3)
 311  FORMAT(1H ,I5,F8.1,1P,22E9.1)
 313  FORMAT(F6.1,1P,33E9.2)
         END
C:::::::::::::::::::::::::::::ECELLS:::::::::::::::::::::::::::::::::::
C...... Subroutine to set energy grid for the 2 stream program
C...... EMAX is the maximum energy required. JMAX=index of max E
C...... It returns the boundaries EB, the midpoint energies E and
C...... the width of the energy cells DELTE
C......  LINEAR energy grid between energies E1 and E2
C...... where the energy steps are S1 and S2
      SUBROUTINE ECELLS(E1,E2,S1,S2,EMAX,JMAX,EB,E,DELTE)
      DIMENSION EB(201),E(201),DELTE(201)
      JMAX=0

      !...... set up cell boundaries EB
      EB(1)=0.0
      DO J=2,199
        IF(EB(J-1).GT.EMAX) GO TO 10
        !... find gradient and abscissa of linear fit for DELTE
        GRAD=(S2-S1)/(E2-E1)
        ABSC=S1-E1*GRAD
        DELTE(J)=(GRAD*EB(J-1)+ABSC)
        IF(DELTE(J).GT.S2) DELTE(J)=S2
        IF(DELTE(J).LT.1.0) THEN
           GRAD=0.0
           ABSC=1.0
           DELTE(J)=1.0
        ENDIF
        
        !.. DELTE must not be less than 1 because the frequencies are 1 eV
        IF(DELTE(J).LT.1.0) DELTE(J)=1.0

        EB(J)=(EB(J-1)+DELTE(J))
        !..WRITE(6,95) J,EB(J),EB(J-1),DELTE(J),(GRAD*EB(J)+ABSC)
      ENDDO

 10   CONTINUE
      JMAX=J-2

      !.. Set up cell midpoint energies E(J) and width DELTE(J)
      DO J=1,JMAX
        E(J)=0.5*(EB(J+1)+EB(J))
        DELTE(J)=EB(J+1)-EB(J)
        !..WRITE(6,95) J,E(J),DELTE(J)
      ENDDO
      EMAX=E(JMAX)

 95   FORMAT(I5,9F9.1)

      RETURN
      END
C::::::::::::::::::::::::::::PITCH:::::::::::::::::::::::::::::::::::
C..... Take pitch angle into account assuming adiabatic variation
C..... along the field line. The pitch angle only varies above the upper
C..... boundary altitude for the PDE calculation Z(M). The program could
C..... be written to avoid this if necessary.
      SUBROUTINE PITCH(J,IS,M,M1,IMAX,T1,T2,PRED,PRODUP,PRODWN,BM,Z)
      DOUBLE PRECISION BM(401),Z(401)
      DIMENSION T1(401),T2(401),PRED(401),PRODUP(201,401)
     > ,PRODWN(201,401)
      !......... RBM used in pitch angle variation
      RBM=0.667/BM(M)
      DO 56 I=M1,IMAX
      !....... pitch angle variation above Z(M)
      IF(Z(I).GT.Z(M)) RAVMU=1.0/SQRT(1.0-BM(I)*RBM)
      IF(Z(I).LT.Z(M)) RAVMU=1.73
      !...... To reset after solution
      IF(IS.EQ.1) RAVMU=1.0/RAVMU
      T1(I)=T1(I)*RAVMU
      T2(I)=T2(I)*RAVMU
      PRED(I)=PRED(I)*RAVMU
      PRODUP(J,I)=PRODUP(J,I)*RAVMU
      PRODWN(J,I)=PRODWN(J,I)*RAVMU
 56   CONTINUE
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C....to calculate cascade production rates from electrons o(1d) and n2*
      SUBROUTINE CASSIM(J,IZ,IWR,ALT,SIGEX,XN,PROCAS,PHIUP,PHIDWN
     > ,E,DELTE,SIGO1D,SIGION,JMAX)
      DIMENSION SIGEX(3),XN(3,401),PROCAS(201,401),PHIUP(401)
     > ,E(201),DELTE(201),PHIDWN(401),SIGION(3)
      DATA BSCX,EBSCX/0.5,0.0/

      !............. cascade from O(1D) and N2*
      PREX=(PHIUP(IZ)*(1-BSCX)+PHIDWN(IZ)*BSCX)*DELTE(J)
      JK=J-2
      IF(JK.GT.0) THEN
         IF(DELTE(J).GT.1) JK=J-1
         CORSIG=2.0/(E(J)-E(JK))
         PROCAS(JK,IZ)=PROCAS(JK,IZ)+PREX*SIGO1D*CORSIG*XN(1,IZ)/
     >   DELTE(JK)
      ENDIF
      IF(E(J).LT.6) PROCAS(J-1,IZ)=PROCAS(J-1,IZ)+PREX*SIGEX(3)*XN(3,IZ)
      !........ O fine structure
      IF(E(J).LE.2) PROCAS(J-1,IZ)=PROCAS(J-1,IZ)+PREX*SIGEX(1)*XN(1,IZ)
        RETURN
        END
C::::::::::::::::::::::::::::CASEXx:::::::::::::::::::::::::::::::::::
C..... to calculate cascade production rates from excitations
C..... JOE= index of bin for depositing degraded electrons
C..... Modified 2003-20-13 to do all 3 species
      SUBROUTINE CASEXx(FLDIM,IE,ISP,IZ,ALT,SIGEX,XN,PROCAS,PHIUP,
     >  PHIDWN,E,DELTE,JOEIN,ELOSS)
      IMPLICIT NONE
      INTEGER FLDIM,IE,ISP,IZ,JOE,JOEM1,JOEM2,JOEIN
      REAL SIGEX(3),XN(3,FLDIM),PROCAS(201,FLDIM),PHIUP(FLDIM)
     > ,E(201),DELTE(201),PHIDWN(FLDIM),SIGION(3)
      REAL BSCX,ALT,PREX,CORSIG,ELOSS
      DATA BSCX/0.5/
      CALL BACKSCAT(E(IE),BSCX)     !.. inelastic backscatter coeff

      JOE=JOEIN
      IF(JOEIN.EQ.IE) JOE=JOEIN-1
      IF(E(IE).LE.ELOSS) RETURN
      JOEM1=JOE
      JOEM2=JOE
      IF(DELTE(JOE).LE.ELOSS.AND.JOE.GE.2) JOEM1=JOE-1
      IF(DELTE(JOE).LE.ELOSS.AND.JOE.GE.2) JOEM2=JOE+1
      IF(JOEM1.LT.1) JOEM1=JOE
      IF(JOEM2.LT.1.OR.JOEM2.GE.IE) JOEM2=JOE
      CORSIG=SIGEX(ISP)*ELOSS/(E(IE)-E(JOE))
      PREX=(PHIUP(IZ)*(1-BSCX)+PHIDWN(IZ)*BSCX)*DELTE(IE)/DELTE(JOE)
      !.. allocate production from bin IE in 2 adjacent bins if small bins
      PREX=0.33333*PREX*CORSIG*XN(ISP,IZ)
      PROCAS(JOE,IZ)=PROCAS(JOE,IZ)+PREX
      PROCAS(JOEM1,IZ)=PROCAS(JOEM1,IZ)+PREX
      PROCAS(JOEM2,IZ)=PROCAS(JOEM2,IZ)+PREX
        RETURN
        END
C::::::::::::::::::::::::::::CASEX:::::::::::::::::::::::::::::::::::
C..... to calculate cascade production rates from excitations
C..... JOE= index of bin for depositing degraded electrons
C..... Modified 2003-20-13 to do all 3 species
      SUBROUTINE CASEX(J,ISP,IZ,IWR,ALT,SIGEX,XN,PROCAS,PHIUP,PHIDWN
     > ,E,DELTE,SIGO1D,SIGION,JMAX,JOEIN,ELOSS)
      DIMENSION SIGEX(3),XN(3,401),PROCAS(201,401),PHIUP(401)
     > ,E(201),DELTE(201),PHIDWN(401),SIGION(3)
      DATA BSCX/0.5/
      CALL BACKSCAT(E(J),BSCX)     !... inelastic backscatter coeff

      JOE=JOEIN
      IF(JOEIN.EQ.J) JOE=JOEIN-1
      IF(E(J).LE.ELOSS) RETURN
      JOEM1=JOE
      JOEM2=JOE
      IF(DELTE(JOE).LE.ELOSS.AND.JOE.GE.2) JOEM1=JOE-1
      IF(DELTE(JOE).LE.ELOSS.AND.JOE.GE.2) JOEM2=JOE+1
      IF(JOEM1.LT.1) JOEM1=JOE
      IF(JOEM2.LT.1.OR.JOEM2.GE.J) JOEM2=JOE
      CORSIG=SIGEX(ISP)*ELOSS/(E(J)-E(JOE))
      PREX=(PHIUP(IZ)*(1-BSCX)+PHIDWN(IZ)*BSCX)*DELTE(J)/DELTE(JOE)
      !.. allocate production from bin J in 2 adjacent bins if small bins
      PREX=0.33333*PREX*CORSIG*XN(ISP,IZ)
      PROCAS(JOE,IZ)=PROCAS(JOE,IZ)+PREX
      PROCAS(JOEM1,IZ)=PROCAS(JOEM1,IZ)+PREX
      PROCAS(JOEM2,IZ)=PROCAS(JOEM2,IZ)+PREX
        RETURN
        END
C:::::::::::::::::::::::::::::::::CASION:::::::::::::::::::::::::::::::::::::
C....... calculate secondary and cascade production from ionizing collisions
      SUBROUTINE CASION(J,IZ,IWR,ALT,XN,PROCAS,PHIUP,PHIDWN
     > ,E,DELTE,ISEC,SECSAV,SIGION,JMAX,JOE,ELOSS,PROB)
      DIMENSION SECSAV(2,401),PROCAS(201,401),XN(3,401),SIGION(3),E(201)
     > ,DELTE(201)
      DATA BSCI/0.5/
      CALL BACKSCAT(E(J),BSCI)     !... inelastic backscatter coeff

      !.. Total production of ions at E(J)
      PRION=SIGION(1)*XN(1,IZ)+SIGION(2)*XN(2,IZ)+SIGION(3)*XN(3,IZ)
      PEFLUP=(PHIUP*(1.0-BSCI)+PHIDWN*BSCI)*DELTE(J)
      !.. total ions produced to energy E(J)
      SECSAV(ISEC,IZ)=SECSAV(ISEC,IZ)+PRION*(PEFLUP)
      PROCAS(J-1,IZ)=PROCAS(J-1,IZ)+SECSAV(ISEC,IZ)*PROB
      !.. degraded primaries

      IF(JOE.LE.1.OR.JOE.GE.JMAX) RETURN
      JOEM1=JOE-1
      JOEP1=JOE+1
      IF(JOEP1.LT.J) THEN
         PRION=0.166666666*PRION*PEFLUP
         PROCAS(JOE,IZ)=PROCAS(JOE,IZ)+2.0*PRION/DELTE(JOE)
         PROCAS(JOEM1,IZ)=PROCAS(JOEM1,IZ)+2.0*PRION/DELTE(JOEM1)
         PROCAS(JOEP1,IZ)=PROCAS(JOEP1,IZ)+2.0*PRION/DELTE(JOEP1)
      ELSE
         PRION=0.3333333*PRION*PEFLUP
         PROCAS(JOE,IZ)=PROCAS(JOE,IZ)+1.5*PRION/DELTE(JOE)
         PROCAS(JOEM1,IZ)=PROCAS(JOEM1,IZ)+1.5*PRION/DELTE(JOEM1)
      ENDIF

      RETURN
      END
C::::::::::::::::::::::::::::::PEPRIM:::::::::::::::::::::::::::::::
C...... this program evaluates  the primary + cascade production rates
C...... by using production frequencies. see richards and torr jgr 1984(5?)
C...... Production frequencies using Samson and Pareek O X-S. Kirby et
C...... al for O2, and N2, and He
      SUBROUTINE PEPRIM(I,XN,HE,PRED,JEMAX,J,EE,DELE,IWR,EFAC)
      INTEGER IDAY
      REAL AP,DEC,ETRAN,BLON,F107,F107A,SEC,ZN2D
      DOUBLE PRECISION COLUM,OTHPR1,OTHPR2,OTHPR3,HE,GL
      DIMENSION XN(3,401),PRED(401)
      !.. Production frequencies from subroutine CONVJS
      COMMON/RJS/RJOX(201),RJN2(201),RJO2(201),RJHE(201)
      COMMON/NPRD/COLUM(3,401),OTHPR1(6,401),OTHPR2(6,401),OTHPR3(6,401)
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)

      EP=EE+17      !..  EP= photon energy
      !.. Set the same primary energy for the 20-22 eV peaks
      IF(EE.GT.20.0.AND.EE.LT.29.0) EP=41.0

      !.. total EUV absorption cross sections are T_XS_??.
      T_XS_O2=2.2*T_XS_OX(EP)    !.. O2 XS is 2.2* O XS

      !.. column densities are from subroutine PRIMPR to get attenuation
      TAU=COLUM(1,I)*T_XS_OX(EP)+COLUM(2,I)*T_XS_O2+
     >  COLUM(3,I)*T_XS_N2(EP)

      AFAC=EXP(-TAU)      !.. EUV attenuation factor
      !.. primary production rates
      PRED(I)=(RJOX(J)*XN(1,I)+RJN2(J)*XN(3,I)+RJO2(J)*XN(2,I)
     > +1.0*RJHE(J)*HE)*1.0E-9*AFAC*EFAC

      RETURN
      END
C::::::::::::::::::::::::::::::::FNDBIN::::::::::::::::::::::::::::
C...... A program to determine which bin to allocate a degraded electron
C...... J= index of the primary of energy EE(J). JMAX is the max J, ELOSS=
C...... energy lost by primary and the index is returned in IE
      SUBROUTINE FNDBIN(J,JMAX,ELOSS,EE,IE)
      INTEGER J,JMAX,IE(JMAX)
      REAL EE(JMAX)
      !.. Energy of degraded primary
      EDPRIM=EE(J)-ELOSS
      IF(EDPRIM.LT.0.0) WRITE(6,*) '   ELOSS too large in FNDBIN'
      IF(EDPRIM.LT.0.0) STOP
      !.. Test degraded energy against next lowest bin until it fits
      DO 29 I=J,2,-1
      IF(EDPRIM.GT.EE(I)) GO TO 30
      IF(ABS(EDPRIM-EE(I)).LE.ABS(0.5*(EE(I)-EE(I-1)))) GO TO 30
 90   FORMAT(I5,F9.2,I5,9F9.1)
 29   CONTINUE
 30   CONTINUE
      !.. store index for return
      IE(J)=I
      IF(ELOSS.LT.0) WRITE(6,'(2I5,9F8.1)') J,IE(J),EE(J),
     >  EE(IE(J)),ELOSS,EE(J)-EE(J-1)
      RETURN
      END
C::::::::::::::::::::::::::::: TRIS1 ::::::::::::::::::::::::::::::
C
C   Modified by Yingtian Long from Dr. Richards' subroutine
C                    /08/27/97
C
C  This subroutine is used to solve photoelectron flux, phi down, by
C  solving the second order PDE. In th calculation for derivatives
C  at a grid point, the variable spaces between adjacent grid points
C  are considered.
C
C... The derivation of the equations was done as a CS 695 course
C... in summer 1997. PGR has a copy of the final paper.
C
      SUBROUTINE TRIS1(IDIR,M2,M,J,MT,BM,BG,Z,IMAX,DH,PRED,PRODUP,
     >  PRODWN)
      DOUBLE PRECISION DH,Z(401),BM(401),BG(401)
      DIMENSION A(401),B(401),C(401),D(401),PRED(401)
     >         ,PRODWN(201,401),PRODUP(201,401) 
 
      COMMON/TRI/PHIDWN(401),PHIUP(401),T1(401),T2(401),DS(401)
      COMMON/PANGS/YPAS,PASK,IPAS,IPASC,FPAS
 
      DATA   ISW, FAC2, JSN, JSW, KITER, M1, ZPRIN, ZZZZ
     >     /  1 ,   0 ,  3 ,   1,   75 ,  1,   300, 1000 /
 
      IMXM1=IMAX-1
      !...... DO LOOP FOR ITERATING THE SOLUTIONS... PHIDWN IS SOLVED USING
      !...... TRIDAG SOLVER FOR ONE HEMISPHERE, PHIUP IS SOLVED ANALYTICALLY TO UPPER
      !...... BDY IN CONJUGATE H-S, PHIUP FOR C.H.S IS FOUND USING TRIDAG, THEN
      !...... PHIDWN IS SOLVED ANALYTICALLY BACK ALONG THE FIELD LINE
      DO 40 I=M2,M
        DELZ=DS(I)+DS(I+1)                            ! ds(i)+ds(i+1)
        DLB=(BM(I+1)-BM(I-1))/(BM(I)*DELZ)            ! (dB)/(B*ds)
        DSLB=2.*((BM(I+1)/DS(I+1)+BM(I-1)/DS(I))/DELZ
     &                  -BM(I)/(DS(I+1)*DS(I)))/BM(I) ! (d^2B)/(B*ds^2)
        DLT1=(T1(I+1)-T1(I-1))/(T1(I)*DELZ)           ! (dT1)/(T1*ds)
        DTS2=(T2(I+1)-T2(I-1))/DELZ                   ! (dT2)/(ds)
        DPR1=(PRED(I+1)-PRED(I-1))/(2*DELZ)           ! 0.5*(dq)/(ds)
        DPR2=(PRODWN(J,I+1)-PRODWN(J,I-1))/DELZ       ! (dq-)/(ds)
        PHI = 1.
        ALPHA = -(DLT1 + 2.*DLB)         ! -[(dT1)/(T1*ds)+2(dB)/(B*ds)]
        BETA = IDIR*(T2(I)*DLT1-DTS2-DSLB)-T2(I)**2+T1(I)**2-ALPHA*DLB
        A(I) = (2.*PHI/DS(I)-ALPHA)/DELZ
        B(I) = BETA-2.*PHI/(DS(I)*DS(I+1))
        C(I) = (2.*PHI/DS(I+1)+ALPHA)/DELZ
        D(I) = (.5*PRED(I)+PRODWN(J,I))*(IDIR*(DLT1+DLB)-T2(I))
     &         -IDIR*(DPR1+DPR2)-T1(I)*(.5*PRED(I)+PRODUP(J,I))
        !!!  where PRED = q, PRODUP = q+ and PRODWN = q-
 40   CONTINUE
 
      !.... END OF D.E. COEFFS --- SOLUTION FOR NEAR H-S ....
      D(M2)=D(M2)-A(M2)*PHIDWN(M2-1)
      D(M)=D(M)-C(M)*PHIDWN(M+1)
      CALL TRIDAG(PHIDWN,M2,M,A,B,C,D)

 
C:::::::: PHIUP IS EVALUATED  ANALYTICALLY   ::::
      DO 60 I=M2,IMXM1
        R1=(T1(I)*PHIDWN(I)+(PRED(I)+2.*PRODUP(J,I))/2.)/T2(I)/BM(I)
        T2DS=T2(I)*DS(I)
        IF(T2DS.GT.70.0)T2DS=70.0

        PHIUP(I)=BM(I)* (R1+(PHIUP(I-1)/BM(I-1)-R1)*EXP(-T2DS))
      !....       IF(J.EQ.31.AND.Z(I).GT.250.AND.Z(I).LT.900)
      !....     > WRITE(55,93) Z(I),T1(I),T2(I),PRED(I),PRODUP(J,I),BM(I),DS(I)
      !....     >  ,R1,PHIDWN(I),PHIUP(I)
        IF(I.EQ.IPAS+1) PHIUP(I)=PHIUP(I)*(1.-FPAS)
 60   CONTINUE
 
      RETURN
 93   FORMAT(1X,F7.0,1P,22E10.1)
      END

C::::::::::::::::::::::::::::::::: TRISM1 :::::::::::::::::::::::::::::::
      SUBROUTINE TRISM1(IDIR,M2,M,J,MT,BM,BG,Z,IMAX,DH,PRED,PRODUP,
     > PRODWN)
      DOUBLE PRECISION DH,Z(401),BM(401),BG(401)
      DIMENSION A(401),B(401),C(401),D(401),PRED(401)
     > ,PRODWN(201,401),PRODUP(201,401)
      COMMON/TRI/PHIDWN(401),PHIUP(401),T1(401),T2(401),DS(401)
      COMMON/PANGS/YPAS,PASK,IPAS,IPASC,FPAS
      DATA   ISW , FAC2 , JSN , JSW , KITER , M1 , ZPRIN , ZZZZ
     >     /  1  ,   0  ,  3  ,  1  ,   75  ,  1 ,  300  , 1000 /
 
      IMXM2=IMAX-MT-1
      !...... DO LOOP FOR ITERATING THE SOLUTIONS... PHIDWN IS SOLVED USING
      !...... TRIDAG SOLVER FOR ONE HEMISPHERE, PHIUP IS SOLVED ANALYTICALLY TO UPPER
      !...... BDY IN CONJUGATE H-S, PHIUP FOR C.H.S IS FOUND USING TRIDAG, THEN
      !...... PHIDWN IS SOLVED ANALYTICALLY BACK ALONG THE FIELD LINE
      DO 40 I=M2,M
        DELZ=DS(I)+DS(I+1)                            ! ds(i)+ds(i+1)
        DLB=(BM(I+1)-BM(I-1))/(BM(I)*DELZ)            ! (dB)/(B*ds)
        DSLB=2.*((BM(I+1)/DS(I+1)+BM(I-1)/DS(I))/DELZ
     &                  -BM(I)/(DS(I+1)*DS(I)))/BM(I) ! (d^2B)/(B*ds^2)
        DLT1=(T1(I+1)-T1(I-1))/(T1(I)*DELZ)           ! (dT1)/(T1*ds)
        DTS2=(T2(I+1)-T2(I-1))/DELZ                   ! (dT2)/(ds)
        DPR1=(PRED(I+1)-PRED(I-1))/(2*DELZ)           ! 0.5*(dq)/(ds)
        DPR2=(PRODUP(J,I+1)-PRODUP(J,I-1))/DELZ       ! (dq+)/(ds)
        PHI = 1.
        ALPHA = -(DLT1 + 2.*DLB)         ! -[(dT1)/(T1*ds)+2(dB)/(B*ds)]
        BETA = IDIR*(T2(I)*DLT1-DTS2-DSLB)-T2(I)**2+T1(I)**2-ALPHA*DLB
        A(I) = (2.*PHI/DS(I)-ALPHA)/DELZ
        B(I) = BETA-2.*PHI/(DS(I)*DS(I+1))
        C(I) = (2.*PHI/DS(I+1)+ALPHA)/DELZ
        D(I) = (.5*PRED(I)+PRODUP(J,I))*(IDIR*(DLT1+DLB)-T2(I))
     &         -IDIR*(DPR1+DPR2)-T1(I)*(.5*PRED(I)+PRODWN(J,I))
        !!!  where PRED = q, PRODUP = q+ and PRODWN = q-
  40  CONTINUE
      !.... END OF D.E. COEFFS --- SOLUTION FOR NEAR H-S ....
      !..*** CONJUGATE SOLUTIONS *******
 10   D(M2)=D(M2)-A(M2)*PHIUP(M2-1)
      D(M)=D(M)-C(M)*PHIUP(M+1)
      CALL  TRIDAG(PHIUP,M2,M,A,B,C,D)
      !.. PHIDWN IS EVALUATED ANALYTICALLY   ::::
      DO 160 I=M1,IMXM2
        K=IMAX-I
        R1=(T1(K)*PHIUP(K)+(PRED(K)+2.*PRODWN(J,K))/2.)/T2(K)/BM(K)
        T2DS=T2(K)*DS(K+1)
        IF(T2DS.GT.70.0)T2DS=70.0
        PHIDWN(K)=BM(K)*(R1+(PHIDWN(K+1)/BM(K+1)-R1)*EXP(-T2DS))
        IF(K.EQ.IPASC-1) PHIDWN(K)=PHIDWN(K)*(1-FPAS)
 160  CONTINUE
 
      RETURN
      END

C:::::::::::::::::::::::::: TRIDAG :::::::::::::::::::::::::::::::::::::
C...... FOR SOLVING A SYSTEM OF LINEAR SIMULTANEOUS EQUATIONS WITH A
C...... TRIDIAGONAL COEFF MATRIX. THE EQNS ARE NUMBERED FROM IF TO L,
C...... & THEIR  SUB-DIAG. , DIAG. ,& SUPER-DIAG COEFFS. ARE STORED
C...... IN THE ARRAYS A, B, C. THE RIGHT HAND SIDE OF THE VECTOR IS
C...... STORED IN D. THE COMPUTED SOLUTION VECTOR IS STORED
C......  IN THE ARRAY DELTA. THIS ROUTINE COMES FROM CARNAHAN, LUTHER,
C......  AND WILKES, APPLIED NUMERICAL METHODS, WILEY, 1969, PAGE 446
      SUBROUTINE TRIDAG(DELTA,IF,L,A,B,C,D)
      DIMENSION A(401),B(401),C(401),D(401)
      DIMENSION ALPHA(401),DELTA(401),GAMMA(401)
      !.....  COMPUTE INTERMEDIATE ARRAYS ALPHA & GAMMA
        ALPHA(IF)=B(IF)
        GAMMA(IF)=D(IF)/ALPHA(IF)
        IFP1=IF+1
        DO 1 I=IFP1,L
        ALPHA(I)=B(I)-A(I)*C(I-1)/ALPHA(I-1)
        GAMMA(I)=(D(I)-A(I)*GAMMA(I-1))/ ALPHA(I)
1       CONTINUE
      !.....  COMPUTE FINAL SOLUTION VECTOR V
        DELTA(L)=GAMMA(L)
        LAST=L-IF
        DO 2 K=1,LAST
        I=L-K
        DELTA(I)=GAMMA(I)-C(I)*DELTA(I+1)/ALPHA(I)
2       CONTINUE
        RETURN
        END
C::::::::::::::::::::: T_XS_N2 :::::::::::::::::::::::::::
C.... This function calculates the N2 total photoionization
C.... cross section. P. Richards 2003-10-04
      REAL FUNCTION T_XS_N2(EP)
      IMPLICIT NONE
      REAL EP   !... photon energy
      REAL ESAVE
      DATA ESAVE/0.0/

      !.. Wavelength < 20 A, Auger ionization
      IF(EP.GE.600.0) THEN              
        T_XS_N2=0.5E-18
      !.. Wavelength < 31 A, Auger ionization
      ELSEIF(EP.GE.400.0) THEN              
        T_XS_N2=1.0E-18
      !.. Wavelength 31.62 to 23.70 A
      ELSEIF(EP.GE.392.0) THEN
        T_XS_N2=EXP(7.9864*ALOG(EP)-91.6604)
      !.. Wavelength 225 to 125 A
      ELSEIF(EP.GE.55.09) THEN
        T_XS_N2=EXP(-2.3711*ALOG(EP)-29.8142) 
      !.. Wavelength > 225 A
      ELSE
        T_XS_N2=EXP(-1.1077*ALOG(EP)-34.8787)  
      ENDIF

      !..IF(NINT(10*EP).NE.NINT(10*ESAVE)) WRITE(6,'(2F8.1,1P,2E10.2)') 
      !..> 12394.224/EP,EP, T_XS_N2/(3.39E-17*EXP(-0.0263*EP)), T_XS_N2
       ESAVE=EP

      !.. old parameterization
      !..T_XS_N2=3.39E-17*EXP(-0.0263*EP)

      RETURN
      END
C::::::::::::::::::::: T_XS_OX :::::::::::::::::::::::::::
C.... This function calculates the OX total photoionization
C.... cross section. P. Richards 2003-10-04
C.... Samson and Pareek Phys. Rev. A, 31, 1470, 1985

      REAL FUNCTION T_XS_OX(EP)
      IMPLICIT NONE
      REAL EP   !... photon energy
      REAL ESAVE
      DATA ESAVE/0.0/

      !.. NEW parameterization
      IF(EP.GE.500.0) THEN                 
        !.. Wavelength shorter than 25 A, Auger ionization
        T_XS_OX=0.5E-18
      ELSEIF(EP.GE.165.26) THEN                 
        !.. Wavelength shorter than 75 A
        T_XS_OX=EXP(-2.5209*ALOG(EP)-28.8855)
      ELSEIF(EP.GE.55.09) THEN              
        !.. Wavelength between 78 and 256.26 A
        T_XS_OX=EXP(-1.7871*ALOG(EP)-32.6335)
      ELSE
        !.. Wavelength longer than 256.26 A
        T_XS_OX=EXP(-1.3077*ALOG(EP)-34.5556)   
      ENDIF

      !..IF(NINT(10*EP).NE.NINT(10*ESAVE)) WRITE(6,'(2F8.1,1P,2E10.2)') 
      !..> 12394.224/EP,EP, T_XS_OX/(27.2E-18*EXP(-3.09E-2*EP)), T_XS_OX
       ESAVE=EP

      !.. old parameterization
      !.. T_XS_OX=27.2E-18*EXP(-3.09E-2*EP)

      RETURN
      END
C::::::::::::::::::::::::::::HE10830XS(EP)::::::::::
C.....Function to calculate the He 10830 cross section
C..... Cross section from Lara Waldrup
      REAL FUNCTION HE10830XS(EP)
      IMPLICIT NONE
      REAL EP   !... photon energy
      REAL pc(0:7)

      !.. fitting coefficients
      data pc/-2.4466E-16,3.2766E-17,-1.7520E-18,5.0229E-20,
     >        -8.4025E-22,8.2195E-24,-4.3589E-26,9.6785E-29/ 
         
      if (EP .lt. 19.82) then
          HE10830XS=0.0
      else if (EP .ge. 19.82 .and. EP .le. 80.0) then
        HE10830XS=pc(0)+pc(1)*EP+pc(2)*EP**2+pc(3)*EP**3+pc(4)*EP**4+
     >            pc(5)*EP**5+pc(6)*EP**6+pc(7)*EP**7
      else
        HE10830XS=7.75E-19*exp(-(EP-100.)/44.)
      endif

      RETURN
      END
